TESS Atlas fit for TOI {{{TOINUMBER}}}¶

Version: {{{VERSIONNUMBER}}}

Note: This notebook was automatically generated as part of the TESS Atlas project. More information can be found on GitHub: github.com/dfm/tess-atlas

In this notebook, we do a quicklook fit for the parameters of the TESS Objects of Interest (TOI) in the system number {{{TOINUMBER}}}. To do this fit, we use the exoplanet library and you can find more information about that project at exoplanet.dfm.io.

From here, you can scroll down and take a look at the fit results, or you can:

Caveats¶

There are many caveats associated with this relatively simple "quicklook" type of analysis that should be kept in mind. Here are some of the main things that come to mind:

The orbits that we fit are constrained to be circular. One major effect of this approximation is that the fit will significantly overestimate the confidence of the impact parameter constraint, so the results for impact parameter shouldn't be taken too seriously.
Transit timing variations, correlated noise, and (probably) your favorite systematics are ignored. Sorry!
This notebook was generated automatically without human intervention. Use at your own risk!

Getting started¶

To get going, we'll need to make out plots show up inline and install a few packages:

%matplotlib inline
!pip install -q -U lightkurve fbpca exoplanet corner pymc3 dynesty isochrones

Then we'll set up the plotting styles and do all of the imports:

import functools
import logging
import multiprocessing as mp
import os
import warnings
from pathlib import Path
from typing import List, Optional
from copy import deepcopy

import arviz as az
from celerite2.theano import terms, GaussianProcess
import corner
import exoplanet as xo
import lightkurve as lk
import matplotlib.pyplot as plt
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import numpy as np
import pandas as pd
import pymc3 as pm
import pymc3_ext as pmx
import theano.tensor as tt
from IPython.display import display
from astroquery.mast import Catalogs

get_ipython().magic('config InlineBackend.figure_format = "retina"')

# TEMPORARY WORKAROUND
try:
    mp.set_start_method("fork")
except RuntimeError: # "Multiprocessing context already set"
    pass
    
# Don't use the schmantzy progress bar
os.environ["EXOPLANET_NO_AUTO_PBAR"] = "true"

# Warning
warnings.filterwarnings("ignore", category=DeprecationWarning)
warnings.filterwarnings("ignore", category=FutureWarning)

# Logging setup
logger = logging.getLogger("theano.gof.compilelock")
logger.setLevel(logging.ERROR)
logger = logging.getLogger("exoplanet")
logger.setLevel(logging.DEBUG)

# matplotlib settings
plt.style.use("default")
plt.rcParams["savefig.dpi"] = 100
plt.rcParams["figure.dpi"] = 100
plt.rcParams["font.size"] = 16
plt.rcParams["font.family"] = "sans-serif"
plt.rcParams["font.sans-serif"] = ["Liberation Sans"]
plt.rcParams["font.cursive"] = ["Liberation Sans"]
plt.rcParams["mathtext.fontset"] = "custom"
plt.rcParams['image.cmap'] = 'inferno'


# Constants
TOI_DATASOURCE = (
    "https://exofop.ipac.caltech.edu/tess/download_toi.php?sort=toi&output=csv"
)

MIN_NUM_DAYS = 0.25

# TOI_NUMBER = {{{TOINUMBER}}}
# __version__ = {{{VERSIONNUMBER}}}
# FILENAME = {{{FILENAME}}}
TOI_NUMBER = 103
__version__ = "TEST"
FILENAME = "template.ipynb"

Fitting stellar parameters¶

Next, we define some code to grab the TOI list from ExoFOP to get the information about the system.

We wrap the information in three objects, a TIC Entry, a Planet Candidate and finally a Lightcurve Data object.

The TIC Entry object holds one or more Planet Candidates (each candidate associated with one TOI id number) and a Lightcurve Data for associated with the candidates. Note that the Lightcurve Data object is initially the same fopr each candidate but may be masked according to the candidate transient's period.
The Planet Candidate holds informaiton on the TOI data collected by SPOC (eg transit period, etc)
The Lightcurve Data holds the lightcurve time and flux data for the planet candidates.

def get_tic_data_from_database(toi_number: int) -> pd.DataFrame:
    """Get rows of about a TIC  from ExoFOP associated with a TOI target.
    :param int toi_number: The TOI number for which the TIC data is obtained
    :return: Dataframe with all TOIs for the TIC which contains TOI {toi_id}
    :rtype: pd.DataFrame
    """
    tois = pd.read_csv(TOI_DATASOURCE)
    toi = tois[tois["TOI"] == toi_number + 0.01].iloc[0]
    tic = toi["TIC ID"]
    tois_for_tic = tois[tois["TIC ID"] == tic].sort_values("TOI")
    if len(tois_for_tic) < 1:
        raise ValueError(f"TOI-{toi_number} data for TIC-{tic} does not exist.")
    return tois_for_tic

class PlanetCandidate:
    """Plant Candidate obtained by TESS."""

    def __init__(self, toi_id: float, period: float, t0: float, depth: float,
                 duration: float):
        """
        :param float toi_id: The toi number X.Y where the Y represents the TOI sub number
        :param float period: Planet candidate orbital period (in days)
        :param float t0: Epoch (timestamp) of the primary transit in Barycentric Julian Date
        :param float depth: Planet candidate transit depth, in parts per million
        :param float duration: Planet candidate transit duration, in days.
        """
        self.toi_id = toi_id
        self.period = period
        self.t0 = t0
        self.depth = depth
        self.duration = duration

    @classmethod
    def from_toi_database_entry(cls, toi_data: dict):
        return cls(
            toi_id=toi_data['TOI'],
            period=toi_data["Period (days)"],
            t0=toi_data["Epoch (BJD)"] - 2457000,  # convert to TBJD
            depth=toi_data["Depth (ppm)"] * 1e-3,  # convert to parts per thousand
            duration=toi_data["Duration (hours)"] / 24.0,  # convert to days
        )
    
    def get_mask(self, t:np.ndarray)->List[bool]:
        """Get mask of when data points in this planet's transit"""
        dur = 0.5 * self.duration
        dur = MIN_NUM_DAYS if dur < MIN_NUM_DAYS else dur
        return np.abs(self.get_timefold(t)) < dur

    def get_timefold(self, t):
        return calculate_time_fold(t, self.t0, self.period)

    def to_dict(self):
        return {
            "TOI": self.toi_id,
            "Period (days)": self.period,
            "Epoch (TBJD)": self.t0,
            "Depth (ppt)": self.depth,
            "Duration (days)": self.duration
        }

def calculate_time_fold(t, t0, p):
    """Function to get time-fold"""
    hp = 0.5*p
    return (t-t0+hp) % p - hp

class LightcurveData:
    """Stores Light Curve data for a single target"""

    def __init__(self, time:np.ndarray, flux:np.ndarray, flux_err:np.ndarray):
        """
        :param np.ndarray time: The time in days.
        :param np.ndarray flux: The relative flux in parts per thousand.
        :param np.ndarray fluex_err: The flux err in parts per thousand.
        """
        self.time = time
        self.flux = flux
        self.flux_err = flux_err
        self.masked = False

    @classmethod
    def from_mast(cls, tic: int):
        """Uses lightkurve to get TESS data for a TIC from MAST"""
        print(f"Searching for lightkurve data with target='TIC {tic}', mission='TESS'")
        search = lk.search_lightcurve(target=f'TIC {tic}', mission="TESS")
        print(f"Downloading {len(search)} observations of lightcurve data (TIC {tic})")
        data = search.download_all()
        print("Completed lightcurve data download")
        data = data.stitch()
        data = data.remove_nans().remove_outliers()
        return cls(
            time=np.ascontiguousarray(data.time.value, dtype=np.float64),
            flux=np.ascontiguousarray(1e3 * (data.flux.value - 1), dtype=np.float64),
            flux_err=np.ascontiguousarray(1e3 * data.flux_err.value, dtype=np.float64)
        )

    def apply_mask(self, transit_mask:List[bool]):
        """Mask lightcurce data to look only at the central "days" duration of data """
        if self.masked:
            raise ValueError("Lightcurve already masked once.")
        len_before = len(self.time)
        self.time = np.ascontiguousarray(self.time[transit_mask])
        self.flux = np.ascontiguousarray(self.flux[transit_mask])
        self.flux_err = np.ascontiguousarray(self.flux_err[transit_mask])
        len_after = len(self.time)
        print(f"Masking reduces lightcurve from {len_before}-->{len_after} points")
        assert len_before >= len_after, f"{len_before}-->{len_after}"
        self.masked = True

class TicEntry:
    """Hold information about a TIC (TESS Input Catalog) entry"""

    def __init__(self, tic: int, candidates: List[PlanetCandidate]):
        self.tic_number = tic
        self.candidates = candidates
        self.lightcurve = None

    @property
    def planet_count(self):
        return len(self.candidates)

    @classmethod
    def generate_tic_from_toi_number(cls, toi: int):
        tois_for_tic_table = get_tic_data_from_database(toi)
        candidates = []
        for index, toi_data in tois_for_tic_table.iterrows():
            candidate = PlanetCandidate.from_toi_database_entry(toi_data.to_dict())
            candidates.append(candidate)
        return cls(
            tic=int(tois_for_tic_table['TIC ID'].iloc[0]),
            candidates=candidates,
        )

    def load_lightcurve(self):
        self.lightcurve = LightcurveData.from_mast(tic=self.tic_number)
        
    def get_combined_mask(self):
        masks = [c.get_mask(self.lightcurve.time) for c in self.candidates]
        return [any(mask) for mask in zip(*masks)]
    
    def mask_lightcurve(self):
        self.lightcurve.apply_mask(self.get_combined_mask())
           
    def to_dataframe(self):
        return pd.DataFrame([candidate.to_dict() for candidate in self.candidates])

    def display(self):
        df = self.to_dataframe()
        df = df.transpose()
        df.columns = df.loc['TOI']
        display(df)

    def setup_outdir(self, version):
        toi = int(self.candidates[0].toi_id)
        output_dir = os.path.join("results", version, toi)
        os.makedirs(output_dir, exist_ok=True)
        self.outdir = output_dir

tic_entry = TicEntry.generate_tic_from_toi_number(toi=TOI_NUMBER)
tic_entry.display()

Now, lets download and plot the TESS light curve data for the Planet Candidates using lightkurve:

tic_entry.load_lightcurve()

Searching for lightkurve data with target='TIC 336732616', mission='TESS'
Downloading 1 observations of lightcurve data (TIC 336732616)
Completed lightcurve data download

def plot_lightcurve_and_masks(tic_entry: TicEntry):
    lc = tic_entry.lightcurve
    fig = make_subplots(rows=2, cols=1, shared_xaxes=True, vertical_spacing=0.02, x_title="Time [days]")
    fig.add_trace(go.Scattergl(
        x=lc.time, y=lc.flux,
        mode='lines+markers', marker_color="black", marker_size=2, line_width=0.1,
        hoverinfo='skip', name='Data'),
        row=1, col=1
    )
    fig.update_yaxes(title_text='Relative Flux [ppt]', row=1,col=1)
    fig.add_trace(go.Scattergl(
        x=lc.time, y=tic_entry.get_combined_mask(),
        mode='lines', line=dict(color="white"), fill='tozeroy', name=f'Combined'),
        row=2, col=1
    )
    for i, candidate in enumerate(tic_entry.candidates):
        fig.add_trace(go.Scattergl(
            x=lc.time, y=candidate.get_mask(lc.time),
            mode='lines',  name=f'Planet {i+1}'),
            row=2, col=1
        ) 

    fig.update_yaxes(title_text='Planet Transiting', row=2,col=1)
    fig.show()

plot_lightcurve_and_masks(tic_entry)

For efficiency purposes, let's extract just the data within 0.25 days of the transits:

tic_entry.mask_lightcurve()

Masking reduces lightcurve from 18084-->2577 points

def plot_masked_lightcurve_flux_vs_time_since_transit(tic_entry:TicEntry, model_lightcurves:Optional[List[float]]=[]):
    num_planets = tic_entry.planet_count
    subplot_titles=[f"Planet {i+1}: TOI-{c.toi_id}" for i, c in enumerate(tic_entry.candidates)]
    fig = make_subplots(rows=num_planets, cols=1, subplot_titles=subplot_titles,vertical_spacing=0.1)
    for i in range(num_planets):
        lc = tic_entry.lightcurve
        planet = tic_entry.candidates[i]
        fig.add_trace(go.Scattergl(
            x=planet.get_timefold(lc.time), y = lc.flux,
            mode='markers',
            marker=dict(
                size=3,
                color=lc.time, 
                showscale=False,
                colorbar=dict(title="Days")
            ),
             name=f"Candidate {i+1} Data"
        ), row=i+1, col=1)
        fig.update_xaxes(title_text="Time [days]", row=i+1, col=1)
        fig.update_yaxes(title_text="Relative Flux [ppt]", row=i+1, col=1)
    for i, model_lightcurve in enumerate(model_lightcurves):
        lc = tic_entry.lightcurve
        planet = tic_entry.candidates[i]
        fig.add_trace(go.Scattergl(
            x=planet.get_timefold(lc.time), y = model_lightcurve,
            mode='markers', name=f"Planet {i+1}"
        ), row=i+1, col=1)
    fig.update_layout(height=300*num_planets)
    fig.update(layout_coloraxis_showscale=False)
    fig.show()

plot_masked_lightcurve_flux_vs_time_since_transit(tic_entry)

That looks a little janky, but it's good enough for now.

The probabilistic model¶

We use the probabilistic model as described in Foreman-Mackey et al 2017 to determine the best parameters to fit the transients present in the lightcurve data.

More explicitly, the stellar light curve $l(t; \vec{\theta})$ is modelled with a Gaussian Process (GP). A GP consists of a mean function $\mu(t;\vec{\theta})$ and a kernel function $k_\alpha(t,t';\vec{\theta})$, where $\vec{\theta}$ is the vector of parameters descibing the lightcurve and $t$ is the time during which the lightcurve is under observation

The parameters describing the lightcurve are $\vec{\theta}$ = {
$p_i$ (orbital periods for each planet),
$d_i$ (transient durations for each planet),
$t0_i$ (transient phase/epoch for each planet),
$b_i$ (impact parameter for each planet),
$r_i$ (planet radius in stellar radius for each planet),
$f0$ (baseline relative flux of the light curve from star),
$u1$ $u2$ (two parameters describing the limb-darkening profile of star)
}

With this we can write $$l(t;\vec{\theta}) \sim \mathcal{GP} (\mu(t;\vec{\theta}), k_\alpha(t,t';\vec{\theta}))\ .$$

Here the mean and kernel functions are:

$\mu(t;\vec{\theta})$: a limb-darkened transit light curve (Kipping 2013)
$k_\alpha(t,t';\vec{\theta}))$: a stochastically-driven, damped harmonic oscillator (SHOTterm)

Now that we have defined our transient model, we can implement it in python:

def build_planet_transient_model(tic_entry):
    n = tic_entry.planet_count
    t0s = np.array([planet.t0 for planet in tic_entry.candidates])
    depths = np.array([planet.depth for planet in tic_entry.candidates])
    periods = np.array([planet.period for planet in tic_entry.candidates])
    
    t = tic_entry.lightcurve.time
    y = tic_entry.lightcurve.flux 
    yerr = tic_entry.lightcurve.flux_err 


    with pm.Model() as my_planet_transient_model:
        ## define planet 𝜃⃗ 
        t0 = pm.Normal("t0", mu=t0s, sd=1.0, shape=n)
        p = pm.Lognormal("p", mu=np.log(periods), sd=1.0, shape=n)
        d = pm.Lognormal("d", mu=np.log(0.1), sigma=10.0, shape=n)
        r = pm.Lognormal("r", mu=0.5 * np.log(depths * 1e-3), sd=1.0, shape=n)
        b = xo.distributions.ImpactParameter("b", ror=r, shape=n)
        planet_parms = [r, d, b]
        
        ## define stellar 𝜃⃗ 
        f0 = pm.Normal("f0", mu=0.0, sd=10.0)
        u = xo.distributions.QuadLimbDark("u")
        stellar_params = [f0, u]
        
        ## define 𝑘(𝑡,𝑡′;𝜃⃗ )  
        sigma = pm.InverseGamma("sigma", alpha=3.0, beta=2 * np.median(yerr))
        rho_gp = pm.Lognormal("rho_gp", mu=0, sd=10)
        sigma_gp = pm.Lognormal("sigma_gp", mu=np.log(np.std(y)), sd=10)
        kernel = terms.SHOTerm(sigma=sigma_gp, rho=rho_gp, Q=1.0 / np.sqrt(2))
        noise_params = [sigma, rho_gp, sigma_gp]

        ## define 𝜇(𝑡;𝜃) (ie light)
        orbit = xo.orbits.KeplerianOrbit(period=p, t0=t0, b=b)
        lightcurves = xo.LimbDarkLightCurve(u).get_light_curve(orbit=orbit, r=r, t=t)
        lightcurve = 1e3 * pm.math.sum(lightcurves, axis=-1) + f0
        lightcurves = pm.Deterministic("lightcurves", lightcurves)
        rho_circ = pm.Deterministic("rho_circ", orbit.rho_star)
        
        # Finally the GP observation model
        residual = y - lightcurve
        gp = GaussianProcess(kernel, t=t, diag=yerr ** 2 + sigma ** 2, mean=lightcurve)
        gp.marginal("obs", observed=y)

        # cache params
        my_params = dict(
            planet_params=planet_parms,
            noise_params=noise_params,
            stellar_params=stellar_params,
        )
    return my_planet_transient_model, my_params, gp


def test_model(model):
    """Test a point in the model and assure no nans"""
    with model:
        test_prob = model.check_test_point()
        test_prob.name="log P(test-point)"
        assert not test_prob.isnull().values.any(), test_prob
        test_pt = pd.Series(
            {k:str(round(np.array(v).flatten()[0],2)) for k,v in model.test_point.items()},
            name="Test Point"
        )
        return pd.concat([test_pt, test_prob], axis=1)

planet_transient_model, params, gp = build_planet_transient_model(tic_entry)
test_model(planet_transient_model)

The test point acts as an example of a point in the parameter space. We can now optimize the model sampling parameters before initialising the sampler.

def get_optimized_init_params(model, planet_params, noise_params, stellar_params):
    """Get params with maximimal log prob for sampling starting point"""
    print("Optimizing sampling starting point")
    with model:
        theta = model.test_point
        kwargs = dict(start=theta, verbose=False, progress=False)
        theta = pmx.optimize(**kwargs, vars=[noise_params[0]])
        theta = pmx.optimize(**kwargs, vars=planet_params)
        theta = pmx.optimize(**kwargs, vars=noise_params)
        theta = pmx.optimize(**kwargs, vars=stellar_params)
    print("Optimization complete!")
    return theta

init_params = get_optimized_init_params(planet_transient_model, **params)

Optimizing sampling starting point
Optimization complete!

Now we can plot our initial model:

def plot_lightcurve_with_inital_model(tic_entry: TicEntry, map_soln):
    lc = tic_entry.lightcurve
    fig = go.Figure()
    make_subplots(rows=2, cols=1, shared_xaxes=True, vertical_spacing=0.02, )
    fig.add_trace(go.Scattergl(
        x=lc.time, y=lc.flux,
        mode='lines+markers', marker_color="black", marker_size=2, line_width=0.1,
        hoverinfo='skip', name='Data'),
    )
    for i in range(tic_entry.planet_count):
        fig.add_trace(go.Scattergl(
            x=lc.time, y=map_soln["lightcurves"][:,i]*1e3,
            mode='lines',  name=f'Planet {i}')
        ) 
    fig.update_layout(
        xaxis_title="Time [days]",
        yaxis_title='Relative Flux [ppt]',
    )
    fig.show()

plot_lightcurve_with_inital_model(tic_entry, init_params)

plot_masked_lightcurve_flux_vs_time_since_transit(
    tic_entry=tic_entry,
    model_lightcurves=[init_params["lightcurves"][:,i]*1e3 for i in range(tic_entry.planet_count)]
)

That looks better!

Now on to sampling:

TUNE = 2000
DRAWS = 2000
CHAINS = 1

np.random.seed(286923464)

def start_model_sampling(model):
    with model:
        trace = pm.sample(
            tune=TUNE,
            draws=DRAWS,
            start=init_params,
            chains=CHAINS,
            cores=1,
            step=xo.get_dense_nuts_step(target_accept=0.9),
        )
        return trace

trace = start_model_sampling(planet_transient_model)

Sequential sampling (1 chains in 1 job)
NUTS: [sigma_gp, rho_gp, sigma, u, f0, b, r, d, p, t0]

Sampling 1 chain for 2_000 tune and 2_000 draw iterations (2_000 + 2_000 draws total) took 154 seconds.
Only one chain was sampled, this makes it impossible to run some convergence checks

Then we can take a look at the summary statistics:

az.from_pymc3(trace)

<xarray.Dataset>
Dimensions:            (b_dim_0: 1, chain: 1, d_dim_0: 1, draw: 2000, lightcurves_dim_0: 2577, lightcurves_dim_1: 1, p_dim_0: 1, r_dim_0: 1, t0_dim_0: 1, u_dim_0: 2)
Coordinates:
  * chain              (chain) int64 0
  * draw               (draw) int64 0 1 2 3 4 5 ... 1995 1996 1997 1998 1999
  * t0_dim_0           (t0_dim_0) int64 0
  * p_dim_0            (p_dim_0) int64 0
  * d_dim_0            (d_dim_0) int64 0
  * r_dim_0            (r_dim_0) int64 0
  * b_dim_0            (b_dim_0) int64 0
  * u_dim_0            (u_dim_0) int64 0 1
  * lightcurves_dim_0  (lightcurves_dim_0) int64 0 1 2 3 ... 2573 2574 2575 2576
  * lightcurves_dim_1  (lightcurves_dim_1) int64 0
Data variables:
    t0                 (chain, draw, t0_dim_0) float64 1.327e+03 ... 1.327e+03
    f0                 (chain, draw) float64 -0.424 -0.2292 ... -0.1887 -0.2112
    p                  (chain, draw, p_dim_0) float64 3.548 3.548 ... 3.548
    d                  (chain, draw, d_dim_0) float64 5.179 0.0002491 ... 1.46
    r                  (chain, draw, r_dim_0) float64 0.09643 ... 0.09751
    b                  (chain, draw, b_dim_0) float64 0.02582 ... 0.01605
    u                  (chain, draw, u_dim_0) float64 0.1475 ... -0.02559
    sigma              (chain, draw) float64 0.3189 0.5328 ... 0.3781 0.44
    rho_gp             (chain, draw) float64 0.05344 0.0386 ... 0.03574 0.04221
    sigma_gp           (chain, draw) float64 0.86 0.9152 1.043 ... 1.012 0.7904
    lightcurves        (chain, draw, lightcurves_dim_0, lightcurves_dim_1) float64 ...
    rho_circ           (chain, draw) float64 1.41 1.41 1.41 ... 1.41 1.41 1.41
Attributes:
    created_at:                 2020-10-12T06:14:30.768914
    arviz_version:              0.9.0
    inference_library:          pymc3
    inference_library_version:  3.9.3
    sampling_time:              154.3821198940277
    tuning_steps:               2000

array([0])

array([   0,    1,    2, ..., 1997, 1998, 1999])

array([0])

array([0])

array([0])

array([0])

array([0])

array([0, 1])

array([   0,    1,    2, ..., 2574, 2575, 2576])

And plot the posterior covariances compared to the values from Pepper et al. (2019):

def plot_posteriors(trace):
    samples = pm.trace_to_dataframe(trace, varnames=["p", "r", "b"])
    corner.corner(samples);

plot_posteriors(trace)

findfont: Font family ['sans-serif'] not found. Falling back to DejaVu Sans.

Finally, we save the posteriors and sampling metadata for future use.

def validate_trace_filename(filename):
    suffix = Path(filename).suffix
    if suffix != ".netcdf":
        raise ValueError(f"{suffix} is an invalid extension.")

def save_trace(trace, filename):
    """Save pymc3 trace as a netcdf file"""
    validate_trace_filename(filename)
    az_trace = az.from_pymc3(trace)
    az_trace.to_netcdf(filename)

def load_trace(filename):
    """Load pymc3 trace from netcdf file and return an arviz InferenceData object"""
    validate_trace_filename(filename)
    return az.from_netcdf(filename)

trace_filename = os.path.basename(FILENAME.replace(".ipynb", ".netcdf"))
save_trace(trace, trace_filename)

Bonus: eccentricity¶

As discussed above, we fit this model assuming a circular orbit which speeds things up for a few reasons. First, setting eccentricity to zero means that the orbital dynamics are much simpler and more computationally efficient, since we don't need to solve Kepler's equation numerically. But this isn't actually the main effect! Instead the bigger issues come from the fact that the degeneracies between eccentricity, arrgument of periasteron, impact parameter, and planet radius are hard for the sampler to handle, causing the sampler's performance to plummet. In this case, by fitting with a circular orbit where duration is one of the parameters, everything is well behaved and the sampler runs faster.

But, in this case, the planet is actually on an eccentric orbit, so that assumption isn't justified. It has been recognized by various researchers over the years (I first learned about this from Bekki Dawson) that, to first order, the eccentricity mainly just changes the transit duration. The key realization is that this can be thought of as a change in the impled density of the star. Therefore, if you fit the transit using stellar density (or duration, in this case) as one of the parameters (note: you must have a different stellar density parameter for each planet if there are more than one), you can use an independent measurement of the stellar density to infer the eccentricity of the orbit after the fact. All the details are described in Dawson & Johnson (2012), but here's how you can do this here using the stellar density listed in the TESS input catalog:

def plot_reweighted_ecentricity_samples(tic_number, trace):
    star = Catalogs.query_object(f"TIC {tic_number}", catalog="TIC", radius=0.001)
    tic_rho_star = float(star["rho"]), float(star["e_rho"])
    print("rho_star = {0} ± {1}".format(*tic_rho_star))

    # Extract the implied density from the fit
    rho_circ = np.repeat(trace["rho_circ"], 100)

    # Sample eccentricity and omega from their priors (the math might
    # be a little more subtle for more informative priors, but I leave
    # that as an exercise for the reader...)
    ecc = np.random.uniform(0, 1, len(rho_circ))
    omega = np.random.uniform(-np.pi, np.pi, len(rho_circ))

    # Compute the "g" parameter from Dawson & Johnson and what true
    # density that implies
    g = (1 + ecc * np.sin(omega)) / np.sqrt(1 - ecc ** 2)
    rho = rho_circ / g ** 3

    # Re-weight these samples to get weighted posterior samples
    log_weights = -0.5 * ((rho - tic_rho_star[0]) / tic_rho_star[1]) ** 2
    weights = np.exp(log_weights - np.max(log_weights))

    # Estimate the expected posterior quantiles
    q = corner.quantile(ecc, [0.16, 0.5, 0.84], weights=weights)
    print("eccentricity = {0:.2f} +{1[1]:.2f} -{1[0]:.2f}".format(q[1], np.diff(q)))

    corner.corner(
        np.vstack((ecc, omega)).T,
        weights=weights,
        plot_datapoints=False,
        labels=["eccentricity", "omega"],
    );

plot_reweighted_ecentricity_samples(tic_entry.tic_number, trace)

rho_star = 0.611063 ± 0.146003
eccentricity = 0.42 +0.27 -0.14

As you can see, this eccentricity estimate is consistent (albeit with large uncertainties) with the value that Pepper et al. (2019) measure using radial velocities and it is definitely clear that this planet is not on a circular orbit.

Citations¶

As described in the :ref:citation tutorial, we can use :func:exoplanet.citations.get_citations_for_model to construct an acknowledgement and BibTeX listing that includes the relevant citations for this model.

with planet_transient_model:
    txt, bib = xo.citations.get_citations_for_model()
print(txt)

This research made use of \textsf{exoplanet} \citep{exoplanet} and its
dependencies \citep{celerite2:foremanmackey17, celerite2:foremanmackey18,
exoplanet:agol20, exoplanet:astropy13, exoplanet:astropy18,
exoplanet:exoplanet, exoplanet:kipping13, exoplanet:luger18, exoplanet:pymc3,
exoplanet:theano}.

print("\n".join(bib.splitlines()[:10]) + "\n...")

@misc{exoplanet:exoplanet,
  author = {Daniel Foreman-Mackey and Rodrigo Luger and Ian Czekala and
            Eric Agol and Adrian Price-Whelan and Timothy D. Brandt and
            Tom Barclay and Luke Bouma},
   title = {exoplanet-dev/exoplanet v0.4.0},
   month = oct,
    year = 2020,
     doi = {10.5281/zenodo.1998447},
     url = {https://doi.org/10.5281/zenodo.1998447}
...

	Test Point	log P(test-point)
t0	1327.25	-0.92
p_log__	1.27	-0.92
d_log__	-2.3	-3.22
r_log__	-2.28	-0.92
b_impact__	-0.19	-1.39
f0	0.0	-3.22
u_quadlimbdark__	0.0	-2.77
sigma_log__	0.21	-0.53
rho_gp_log__	0.0	-3.22
sigma_gp_log__	1.57	-3.22
obs	NaN	-6473.11

TOI	103.01
TOI	103.010000
Period (days)	3.547854
Epoch (TBJD)	1327.252563
Depth (ppt)	10.424372
Duration (days)	0.145597